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Microscopic input to a universal nuclear energy density functional can be provided through the 
density matrix expansion (DME), which has recently been revived and improved. Several DME 
implementation strategies are tested for neutron drop systems in harmonic traps by comparing 
to Hartree-Fock (HF) and ab initio no-core full configuration (NCFC) calculations with a model 
interaction (Minnesota potential). The new DME with exact treatment of Hartree contributions is 
found to best reproduce HF results and supplementing the functional with fit Skyrme-like contact 
terms shows systematic improvement toward the full NCFC results. 



PACS numbers: 21.10.-k,21.30.-x,21.60.Jz 

I. INTRODUCTION 

Experiments at radioactive ion beam facilities, stud- 
ies of astrophysical systems such as neutron stars and 
supernovae, and nuclear energy and security needs have 
motivated multipronged efforts to develop nuclear energy 
density functionals (EDF's) with substantially reduced 
errors and improved predictive power away from stabil- 
ity. While great progress has been made in extending 
the reach of ab initio wave function methods beyond the 
lightest nuclei 

[lHl, the EDF approach remains the most 
computationally feasible method for a comprehensive de- 
scription of medium and heavy nuclei @. However, the 
ab initio methods are vital tools for reaching the goal of 
robust functionals informed by microscopic internucleon 
interactions. As part of an ongoing program to achieve 
this goal, in this paper we investigate trapped neutron 
drops with a model interaction. In particular, EDF cal- 
culations using several density matrix expansion (DME) 
implementations are confronted with ab initio no-core full 
configuration (NCFC) 0] results. 

The comparisons presented here exploit developments 
achieved within the Universal Nuclear Energy Density 
Functional (UNEDF) SciDAC-2 collaboration ]3, Q . The 
UNEDF project aims to develop a comprehensive the- 
ory of nuclear structure and reactions utilizing the most 
advanced computational resources and algorithms avail- 
able, including high-performance computing techniques 
to scale to petaflop platforms and beyond Q. One ele- 
ment of the UNEDF program involves the direct injec- 
tion of microscopic physics into novel energy functionals, 
with the DME a key tool @-[l2[. Another element has 
led to efficient density functional theory (DFT) solvers 
adapted for the DME [lj| and to neutrons in external 
traps, which allow accurate and rapid testing of candi- 
date functionals [lij] . A third element is the extensive de- 
velopment of ab initio methods, including improved com- 



putational efficiencies [15[ and extrapolation techniques 
for the NCFC that allow exact calculations (with er- 
ror bars) of the same neutron drop systems to which the 
functionals are applied. 

Neutron drops are a powerful theoretical laboratory 
for improving existing nuclear energy functionals. Micro- 
scopic input to EDF's is particularly needed for neutron- 
rich nuclei, where there are fewer constraints from exper- 
iment and where novel density dependencies from long- 
range pion exchange may be particularly important to 
incorporate. The properties of uniform neutron matter 
have been used in the past as a constraint on phenomeno- 
logical functionals (e.g., see Refs. [jl [13), but compu- 
tational advances now allow accurate microscopic many- 
body calculations of inhomogeneous neutron drops in ex- 
ternal potentials using quantum Monte Carlo or NCFC 
methods [TH, [l9[. These calculations can be used to 
identify deficiencies in existing functionals (e.g., of the 
Skyrme type as in Ref. [TH), to suggest or calibrate new 
versions, or simply to provide control data that supple- 
ments experiment for the optimization of functionals. 

Neutron drops also provide favorable environments for 
the development and testing of non-empirical (i.e., mi- 
croscopically based) functionals. The necessity of an ex- 
ternal potential (because the untrapped system is un- 
bound, with positive pressure) is turned into a virtue by 
allowing external control over the environment. Density 
Functional Theory (DFT) , which provides the theoretical 
underpinning and computational framework for building 
a nuclear EDF, dictates that the same functional applies 
for any external potential, which can therefore be varied 
to probe and isolate different aspects of the EDF. In con- 
trast, the treatment of self-bound systems (such as ordi- 
nary nuclei) has much less flexibility. Furthermore, there 
are serious complications from symmetry breaking [2fj| . 
particularly for the relatively small systems where ab ini- 
tio methods can also be applied. 



2 



The density matrix expansion was introduced long ago 
by Negele and Vautherin, who applied it to G-matrix ef- 
fective interactions to derive a Skyrme-likc Hartree-Fock 
(HF) energy density functional. The DME has been re- 
visited in recent years with the goal of applications to nu- 
clear interactions sufficiently soft that many-body pertur- 
bation theory (MBPT) for nuclei is a quantitative frame- 
work. There are new formal DME developments, as well 
as new formulations that include hybrids between purely 
ab initio and phenomenological functionals. These must 
be tested and validated (or discarded if found to be inad- 
equate); in general, due to the complexity of the various 
DME procedures and the involved density dependence 
they generate for the subsequent functional, we have no 
a priori guidance for the accuracy. 

In this paper, we isolate the DME issues by using a 
simplified model interaction, the two-body Minnesota po- 
tential. Focusing on neutron drops is also advantageous 
because the self-consistent solution of a self-bound sys- 
tem magnifies approximation errors to the extent that it 
is difficult to analyze them. Of course, the main draw- 
back of using neutron drops is the absence of control by 
experimental data. In our case this is not so relevant 
as we use a simple model interaction sufficient for con- 
ducting our basic tests of these methods. Our results 
here provide a foundation for testing more realistic in- 
teractions and improved functionals (e.g., that include 
pairing), and for extensions to self-bound nuclei. 

In Section [Hj we briefly review the methods and in- 
puts used. Results for the Minnesota potential are given 
in Section IIIII for representative trap potentials and two 
closed-shell neutron drop systems. We summarize our 
observations and discuss the next steps to take in Sec- 
tion El 



II. BACKGROUND 

The present calculations combine ingredients from sev- 
eral parts of the UNEDF project. The technical details 
are described in full elsewhere, so we merely review the 
essential features. 



A. DME 

The density matrix expansion (DME) introduced by 
Negele and Vautherin [2f|, [22j provides a route to an 
EDF based on microscopic nuclear interactions through 
a quasi-local expansion of the energy in terms of var- 
ious densities: p(R), t(R), V 2 ,o(R), and so on. Kohn- 
Sham single-particle potentials are immediately obtained 
from simple functional derivatives with respect to these 
densities. The DME originated as an expansion of 
the Brueckner- Hartree-Fock en ergy constructed using the 
nucleon-nucleon G-matrix [U [22[ , which was treated in 
a local (i.e., diagonal in coordinate representation) ap- 
proximation. 



The DME has been reformulated for spin-saturated nu- 
clei using non-local low-momentum interactions in mo- 
mentum representation Q, for which G-matrix summa- 
tions are not needed because of the softening of the in- 
teraction. When applied to a Hartree-Fock energy func- 
tional, the DME yields an EDF in the form of a gener- 
alized Skyrmc functional that is compatible with exist- 
ing codes, by replacing Skyrme coefficients with density- 
dependent functions. As in the original application, a 
key feature of the DME is that it is not a pure short- 
distance expansion but includes resummations that treat 
long-range interactions correctly in a uniform system. 

Extensions of the first calculations from Ref. have 
modified the DME formalism from Negele and Vau- 
therin [2f| , which provides an extremely poor description 
of the vector part of the density matrix [ll| . The stan- 
dard DME is much better at reproducing the scalar den- 
sity matrices, but errors arc still sufficiently large that 
the discrepancy with full finite-range Hartree-Fock cal- 
culations can reach the Me V p er particle level. Gebre- 
mariam and collaborators jllj introduced a new phase 



space averaging (PSA) approach that accounts for the 
diffuse Fermi surface [23| and anisotropy 24[ of the local 
momentum distribution, with no free parameters. The 
PSA treatment leads to substantial improvements, par- 
ticularly for the vector density matrices, where relative 
errors in integrated quantities are reduced by as much as 
an order of magnitude across isotope chains flll |. I n the 
present work, we test the difference between the original 
Negele- Vautherin (NV) and the new PSA prescriptions 
only for scalar parts. 

The Fock energy exhibits spatial non-localities due to 
the convolution of finite-range interaction vertices with 
non-local density matrices. The DME factorizes the non- 
locality of the one-body density matrix by expanding it 
into a finite sum of terms that are separable in rela- 
tive r = T\ — r-i and center of mass R = (ri + r*2)/2 
coordinates. For example, in notation introduced in 
Refs. [Hldli], one expands the spin-scalar part (in both 
isospin channels) of the one-body density matrix as 



p t {n,r 2 ) « J2 n„(fcr) V n (R) , 



(1) 



n=0 



where the functions {V n (R)} denote various local den- 
sities and their gradients (through second order in 
the present work) and H n (kr) denotes the so-called 
n— functions, which depend on the particular formula- 
tion of the DME (here NV or PSA). The arbitrary mo- 
mentum k sets the scale for the decay in the off-diagonal 
direction. As in Ref. we will take k to be the local 
Fermi momentum related to the isoscalar density through 



k F (R) 



3tt s 



po(R] 



1/3 



(2) 



although other choices are possible that include addi- 
tional r and Ap dependencies (25| . 
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It is possible to apply the DME to both Hartrce and 
Fock energies so that the complete Hartree-Fock energy 
is mapped into a local functional. From the earliest 
DME work, however, it was found that treating the 
Hartree contributions exactly provides a better reproduc- 
tion of the density fluctuations and the energy produced 
from an exact HF calculation [22|, [2(|. Restricting the 
DME to the exchange contribution significantly reduces 
the self-consistent propagation of errors [22j . In addi- 
tion, treating the Hartrce contribution exactly does not 
complicate the numerical solutions of the resulting self- 
consistent equations compared to applying the DME to 
both Hartree and Fock terms. Here we will compare sev- 
eral prescriptions for handling the Hartree contribution, 
including a Taylor series expansion (27j . 

A consistent and systematic extension of the DME 
procedure beyond the Hartree-Fock level of MBPT has 
yet to be formulated. For now, attempts to micro- 
scopically construct a quantitative Skyrme-like EDF use 
some ad hoc approximations (e.g., using averaged rather 
than state-dependent energy denominators) when apply- 
ing the DME to iterated contributions beyond the HF 
level and/or re- introduce some phenomenological param- 
eters to be adjusted to data [U [22|, [28l - [3G |. Recent 
work based on chiral NN and NNN effective field the- 
ory interactions motivates an approach to building upon 
the DME/HF functional that incorporates important mi- 
croscopic physics while exploiting the highly developed 
Skyrme EDF technology [nl-GJ]. 

The structure of the chiral interactions is such that 
each coupling in the DME/HF functional is decom- 
posed into a density-independent coupling constant aris- 
ing from zero-range contact interactions and a coupling 
function of the density arising from the universal long- 
range pion exchanges. This clean separation between 
long- and short-distance physics at the HF level sug- 
gests a semi-phcnomcnological approach where the long- 
distance couplings (g™ (R; V-n)) are kept as is, and the 
zero-range constants C t m are optimized to finite nuclei 
and infinite nuclear matter properties [Til fl2j. Thus, 

ST = 9 r(R;V v ) + Cr(V ct ), (3) 

and so on, so that the DME functional splits into two 
terms, 

E[p}=E ct [p}+E 7I [p], (4) 

where the first term E ct [p] collects all contributions from 
the contact part of the interaction plus higher-order 
short-range contributions encoded through the optimiza- 
tion to nuclei and nuclear matter, while the second term 
En [p] collects the long-range NN and NNN pion exchange 
contributions at the Hartree-Fock level. 

Because the contact contributions have essentially the 
same structure as those entering empirical Skyrme func- 
tionals, a microscopically guided Skyrme phenomenology 
has been suggested in which the contact terms in the 
DME functional are released for optimization to finite- 
density observables (Til [l2j. This empirical procedure 



is supported by the observation that the dominant bulk 
correlations in nuclei and nuclear matter are primar- 
ily short-ranged in nature, as evidenced by Brueckner- 
Hartree-Fock (BHF) calculations where the Brueckner 
G-matrix "heals" to the free-space interaction at suffi- 
ciently large distances. One can loosely interpret the re- 
fit of the Skyrme constants to data as approximating the 
short-distance part of the G-matrix with a zero-range 
expansion through second order in gradients. In doing 
so, the constants can capture short-range correlation en- 
ergy contributions beyond Hartree-Fock. We will test 
this strategy for incorporating BHF correlations from the 
Minnesota potential, with the free parameters of the vol- 
ume part of the functional fixed to properties of infinite 
neutron matter and the free surface parameter adjusted 
to NCFC results. We will also consider a direct density- 
dependent modification to model BHF correlations. 



B. Minnesota Potential 

All of the calculations reported here use the Minnesota 
potential. This is a local NN-only potential that is the 
sum of three Gaussians in the radial coordinate (3lj : 

Vij = [VR + lil + P^Vt + ^l-P^V^il + P^ , (5) 

where P a and P r are spin and space exchange operators, 
respectively, and 

Vr = V 0R e~ KRr % , (6) 
V t = -V ot e- K ^ , (7) 
V s = -V 0s e-^ . (8) 

The parameters defining the Vij are given in Table |U 
(Note: we have taken the exchange-mixture parameter u 
in Ref. [3l[ equal to be one.) The Minnesota potential 
reproduces NN effective range parameters and gives rea- 
sonable results for the binding energies of light nuclei. It 
is often used as a semi-realistic potential in model calcu- 
lations. 

TABLE I: Parameters defining the Minnesota potential, see 
Eqs. ©-©. 





value 


Ka 


value 


V 0R 


200.0 MeV 


KR 


1.487 mi -2 


V ot 


178.0 MeV 


Kl 


0.639 fm~ 2 


v 0s 


91.85 MeV 


K>s 


0.465 fm~ 2 



For our purposes, the important characteristics of this 
potential are that it is local, which makes possible the 
immediate adaptation of current DME technology (which 
is not fully developed for non-local potentials), and that 
it is moderately soft, so that HF is a reasonable starting 
point and convergence is adequate in the NCFC. Because 
we plan to use low-momentum interactions in the future, 
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this softness is consistent rather than a shortcoming of 
the model. While the Minnesota potential lacks impor- 
tant features of realistic interactions, such as tensor forces 
and three-nucleon interactions, it provides a convenient, 
non-trivial test case for the DME that sets the stage for 
future tests. 



C. EDF Solvers 

The DME-bascd functionals described in the last sec- 
tion and in Ref . [13 1 have been implemented in the DFT 
solvers HFBRAD|3l and HFBTHO [H- HFBRAD 
is a very fast solver for spherical nuclei and density- 
dependent local density approximations, while HFBTHO 
is much slower but calculates spherical and axially de- 
formed nuclei, and can handle additional gradient cor- 
rections. A Fortran module for both solvers has been 
developed to implement the density-dependent parts of 
the EDF from the DME applied to chiral effective po- 
tentials [34]] . The module contains all of the lengthy 
expressions for the DME couplings and their functional 
derivatives with respect to the density matrix, and for 
numerically stable approximations. The module also has 
the capability to calculate related infinite nuclear matter 
properties. We have developed a similar module that can 
handle expressions coming from the DME of the Min- 
nesota potential plus external potentials. This module 
was linked to existing DFT solvers to calculate results 
presented here. 



D. NCFC 

To test the DME calculations, we use the ab initio no- 
core full configuration (NCFC) approach [(J HU to pro- 
vide exact results (with errors bars) . The NCFC is closely 
related to the no-core shell model (NCSM) HHH, as 
both employ a many-body harmonic oscillator basis that 
treats all nucleons as spectroscopically active. The ba- 
sis space includes all many-body states with excitation 
quanta less than or equal to A r max . A general feature is 
the possibility to completely remove spurious center-of- 
mass excitations, but the present application with an ex- 
ternal potential does not exploit this capability. Rather, 
the center-of-mass motion is part of the physical system. 
The main difference between the NCFC and NCSM is 
that the NCSM employs an interaction renormalized to 
the finite many-body basis, such as the Lee-Suzuki effec- 
tive interaction. 

The NCFC approach involves the extrapolation of a se- 
quence of finite matrix results with the bare interaction 
(as opposed to a Lee-Suzuki effective interaction) to the 
infinite basis space limit. This makes it possible to ob- 
tain basis-space-independent results for binding energies 
and other obscrvables, and to evaluate their numerical 
uncertainties. The extrapolation methods are described 
in Ref. @. A recent calculation of 14 F in Ref. Q, made 



prior to the first experimental measurements, illustrates 
the predictive power of the NCFC approach when cou- 
pled with leadership class computer resources. Note that 
the present calculations do not exploit the full capabili- 
ties of the codes and computers available to further min- 
imize theoretical errors, because current error bars are 
small enough for the present application. 

To solve for ground-state energies, radii, and form fac- 
tors of trapped neutron systems, the code MFDn poMi^ 
was generalized to allow for external potentials. Only 
spherically symmetric harmonic oscillator traps are used 
in the current investigation, but other shapes and de- 
formed traps are also directly available. 



III. RESULTS 

The neutrons are confined by an external single- 
particle harmonic potential: 



(9) 



with harmonic oscillator parameter Ml varied from 
5 MeV to 20 MeV. The calculations here use N = 8 and 
N = 20 neutrons, which form closed shells. In the fu- 
ture we will revisit this problem with pairing included 
and consider intermediate N values. Accurate NCFC re- 
sults are limited to larger oscillator parameters because 
of slow convergence with iV max for the Minnesota po- 
tential. (Note: quantum Monte Carlo techniques such as 
GFMC or AFDMC are effective for smaller Ml and could 
be used for additional comparisons.) In some cases, ex- 
trapolations are not reliable and so only upper bounds 
to the total energy are given. 

Comparisons between different DME treatments of the 
Hartree term are given in Tables [H] IIIII and IIVI The 
full HF results provide a baseline for comparison of the 
Negele-Vautherin (HF/NV) and Phase-Space- Averaging 
(HF/PSA) approximations to HF, with variations based 
on how the Hartree part of the DME is treated. For 
each of the two DME implementations, there are three 
possibilities: treat Hartree with the same DME (NV or 
PSA), use a naive Taylor expansion (NT), or treat it 
exactly. We split the total energy into internal and trap 
contributions, with 



in I 



Etnt ~ U c - 



and 



47r / drr 2 v cyit (r)p(r) 



(10) 



(11) 



Results are presented as deviations from the full HF re- 
sults of the total and internal energies and the radii. 
These results are for spherical solutions, which were 
shown to minimize the energy. That is, by imposing a 
non-zero quadrupolc moment as a constraint, we found 
in all cases that the total energy rapidly increases as the 
quadrupolc moment deviates from zero. 
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TABLE II: Comparison of DME approximations to HF to- 
tal energies for different treatments of the Hartree term, 
expressed as deviations from the full Hartree-Fock results 
AE, = Ei- Ehf in MeV. 







HF/NV 




HF/PSA 




TV 


m 


NV NT exact 


PSA 


NT exact 


8 


3 


0.1 0.2 


0.1 


0.0 


0.1 


0.0 


8 


5 


0.4 0.8 


0.4 


-0.1 


0.0 


0.2 


<s 


10 


2.1 5.1 


2.0 


-1.7 


4.1 


0.9 


8 


15 


4.2 12.9 


4.6 


-7.1 


10.8 


2.1 


8 


20 


6.0 24.2 


7.7 




20.9 


3.4 


20 


3 


0.5 0.8 


0.6 


-0.1 


0.4 


0.2 


20 


5 


1.8 3.4 


2.3 


-1.0 


2.0 


0.9 


20 


10 


5.9 18.5 


11.0 


-14.0 


12.0 


3.9 


20 


15 


3.8 44.3 


22.7 




31.6 


7.9 


20 


20 


-17.8 80.0 


34.8 




61.3 


12.5 



TABLE III: Comparison of DME approximations to HF in- 
ternal energies for different treatments of the Hartree term, 
expressed as deviations from the full Hartree-Fock results 
AE, = Ei- Ehf in MeV. 







HF/NV 




HF/PSA 




TV 


m 


NV NT exact 


PSA NT exact 


8 


3 


0.0 -0.0 


0.0 


0.1 0.0 


0.1 


8 


5 


-0.1 -0.4 


-0.1 


0.2 -0.3 - 


0.1 


8 


10 


0.1 -1.0 


-0.2 


1.3 -1.0 - 


-0.1 


8 


15 


1.1 -1.5 


0.3 


5.4 -1.8 


0.1 


8 


20 


3.2 -1.9 


1.1 


-2.7 


0.5 


20 


3 


-0.2 -0.4 - 


-0.3 


0.1 -0.3 - 


-0.1 


20 


5 


-0.3 -1.2 


-0.6 


1.0 -0.8 - 


-0.2 


20 


10 


3.3 -2.2 


0.3 


11.9 -2.6 


0.2 


20 


15 


16.9 -2.0 


4.4 


-5.4 


1.4 


20 


20 


68.3 -1.2 


11.1 


-8.7 


3.2 



It is evident that the DME with PSA and exact Hartree 
is systematically the closest to HF energies and radii. For 
all DME approximations, internal energies are generally 
closer to full HF than the total energies. This can be un- 
derstood because the internal energy comes mostly from 
the volume part of the EDF. DME and full HF agree 
for infinite (uniform) matter, so similar results can be 
expected from this part. However, one can see that the 
DME approximations give slightly larger radii, which im- 
plies slightly larger external energies, and so larger total 
energies. A possible explanation is that the DME cou- 
pling for the pAp term may not take into account all the 
surface effects. 

Some of the DME/PSA entries in Tables HQ and [IV] 
are missing. In these cases, the calculation failed to con- 
verge because of EDF instabilities [44[ . The instabilities 
occurred at high values of Ml when both the Hartree and 
Fock terms were taken from the DME. This could indi- 
cate some problems at high density in the Hartree part of 
the DME expressions. These instabilities are not, how- 
ever, just simply related to the infinite neutron matter 
properties, and are therefore more involved |44j |. 

An alternative to the DME for microscopically based 
EDF's uses the more completely microscopic but com- 



TABLE IV: Comparison of DME approximations to HF rms 
radii for different treatments of the Hartree term, expressed 
as deviations from the full Hartree-Fock results Art — n — thf 
in fm. 







HF/NV 




HF/PSA 


TV 


m 


NV 


NT exact 


PSA NT 


exact 


8 


3 


0.01 


0.02 


0.00 


-0.02 0.01 


-0.01 


8 


5 


0.03 


0.07 


0.03 


-0.01 0.05 


0.01 


8 


10 


0.04 


0.11 


0.04 


-0.06 0.09 


0.02 


8 


15 


0.03 


0.14 


0.04 


-0.13 0.12 


0.02 


8 


20 


0.02 


0.16 


0.05 


0.15 


0.02 


20 


3 


0.02 


0.05 


0.03 


-0.02 0.02 


0.01 


20 


5 


0.04 


0.09 


0.06 


-0.04 0.06 


0.02 


20 


10 


0.02 


0.13 


0.07 


-0.18 0.09 


0.02 


20 


15 


-0.04 


0.16 


0.07 


0.13 


0.03 


20 


20 


-0.20 


0.18 


0.05 


0.15 


0.02 



putationally far more intensive Optimized Effective Po- 
tential (OEP) method [T3|. In Ref. [Hj], the OEP was 
applied to the same model problem of the Minnesota po- 
tential for neutrons in a trap. Comparisons made to exact 
HF results show that the exact exchange version of OEP 
is almost indistinguishable from HF, in contrast to the 
small but non-negligible discrepancies found here. Fu- 
ture comparisons as both methods continue to be refined 
will help to gauge the accuracy of DME approximations 
and guide the development of corrections. 

To test schemes for incorporating correlations beyond 
HF, we use BHF calculations of neutron matter, which we 
expect to be quite accurate for the Minnesota potential. 
Two strategics are considered, following the discussion in 
Section |HAl 

The first strategy is based on the empirical observa- 
tion that the ratio of the neutron matter HF and BHF 
results is a rather smooth function of density, which we 
denote /(fcp). By assuming a rank-one separable expan- 
sion of the potential, V(k, k'), the G-matrix would take 
the form G(k, k') ~ V(k, \c')/f(kp) and then taking a 
simple Gaussian for the potential form factor and ex- 
panding out the integral that appears in the definition of 
/(fcp), one motivates the form 

/ = a + bp 1/3 + cp+ dp 5 / 3 + ... . (12) 

The coefficient a, 6, c, and d are determined by a fit. 
(Calculations omitting d were also made and yield sim- 
ilar results except for the densest neutron drops.) This 
strategy is implemented in the DFT solvers by evaluat- 
ing the p dependence in the Fock terms as p — > p(R), 
which means the DME/HF couplings simply get scaled 
by f(p). In the exact Hartree treatment, the prescrip- 
tion p —> l/2[p{r\) + pfo)] is used and otherwise p(R) is 
used. This approach is labeled B HF/PSA (or just BHF) 
in the subsequent tables and figures (only results based 
on exact-Hartree, DME/PSA are given). 

The second strategy follows Ref. [13[ to incorporate 
BHF correlations by adding a contact part to the HF 
functional for the Minnesota potential. In general, the 
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contact part E ct [p] of the EDF has the form of the stan- 
dard Skyrme functional 

H c t(r) = ^T +U$(r)+W?(r), (13) 

where 

H?(r) = {C& + cf D pl)pi + Crptn + Cf/H^Pt 
+ C t pvJ p t VJ t + C t j2 J? , (14) 

and the isospin index t = {0, 1} labels isoscalar and 
isovector densities, respectively. In analogy to Ref. [Hj], 
the neutron coupling constants C 4 = Cq + C[, i = 
{p 2 , pr} are fitted to reproduce the neutron matter BHF 
results. This allows us to constrain the zero-range volume 
parameters of the DME-based functional, but not the pa- 
rameters entering the surface part of the functional. 

We optimize the surface parameters in the DME func- 
tional in a manner similar to the optimization done for 
standard Skyrme functionals. One could think of pro- 
cedures based on semi-infinite neutron properties. or 

on the leptodcrmous expansion of the functional [46|. 
Here we optimize the surface coupling constant C pAp = 
C£ Ap + C pAp to NCFC E tot values presented in Tables 
IVl and IVII by using theoretical error bars as weights. A 
simple minimization of the root-mean-square deviation 
yields almost the same result. The values for the neu- 
tron parameters of the Minnesota model are: 

C p2 = -18.25 MeVfm 3 , C pT = 4.57McVfm 5 , 
C pAp = -1.8MeVfm 5 , (15) 

with = C pVJ = C j2 = 0. Calculations done with 
these parameters are labeled as "fit/PSA" in Tables IVl 
and IVII and "fit" in the figures. 

2 

We have set the coupling constant C D to zero in the 
fit to neutron matter properties. In the usual Skyrme- 
DFT scheme, the density dependence controlled by the 
power 7 is needed to produce reasonable saturation prop- 
erties. For example, the incompressibility of symmetric 
nuclear matter is strongly affected by 7, and usually an 
acceptable value requires 7 < 1. In the present calcu- 
lations we do not constrain symmetric matter; indeed, 
the Minnesota potential does not produce realistic sat- 
uration. More generally, the density dependence from 
2 

nonzero C D is used to effectively account for beyond-HF 
and three-body effects. The simplicity of the NN-only 
Minnesota potential seemingly lets us transfer beyond- 
HF effects to the other coupling constants and omit the 

2 

C D term entirely. 

In Tables [V] and IVII we summarize results for different 
trap parameters from the NV and PSA DME implemen- 
tations of HF with exact Hartree, the two strategies to 
go beyond HF (BHF and fit, both based on PSA with ex- 
act Hartree), along with full HF and NCFC results. As 
already observed from the earlier tables, the comparisons 



TABLE V: Results for calculations of 8 neutron drops in har- 
monic potentials. All energies are in MeV and the rms radii 
r lms are in fm. The NCFC results use up to the iV m „ in 
square brackets and parenthesis indicate the extrapolation un- 
certainty in the last quoted digit (s). The approximations are 
explained in the text. 



Approx. 


m 


Etot 


-Bint 






HF 


5 


71.9 


37.6 


34.3 


3.78 


HF/NV 


5 


72.3 


37.5 


34.8 


3.80 


HF/PSA 


5 


72.1 


37.5 


34.5 


3.78 


BHF/PSA 


5 


68.8 


34.9 


33.9 


3.75 


fit/PSA 


5 


70.0 


36.8 


33.2 


3.71 


NCFC [14] 


5 


< 69.5 








HF 


10 


142.4 


69.6 


72.8 


2.75 


HF/NV 


10 


144.5 


69.5 


75.0 


2.79 


HF/PSA 


10 


143.4 


69.6 


73.8 


2.77 


BHF/PSA 


10 


139.4 


66.2 


73.2 


2.75 


fit/PSA 


10 


138.6 


67.3 


71.3 


2.72 


NCFC [16] 


10 138.1(6) 


66(2) 


72(2) 


2.73(3) 


HF 


15 


217.4 


101.8 


115.6 


2.31 


HF/NV 


15 


222.1 


102.1 


120.0 


2.35 


HF/PSA 


15 


219.5 


101.9 


117.6 


2.33 


BHF/PSA 


15 


214.8 


98.5 


116.2 


2.31 


fit/PSA 


15 


212.5 


98.1 


114.4 


2.30 


NCFC [16] 


15 212.7(2) 


98.6(4) 


114.1(4) 2.293(4) 


HF 


20 


296.4 


135.1 


161.3 


2.04 


HF/NV 


20 


304.1 


136.3 


167.8 


2.09 


HF/PSA 


20 


299.8 


135.6 


164.2 


2.06 


BHF/PSA 20 


294.1 


131.8 


162.4 


2.05 


fit/PSA 


20 


290.9 


130.0 


160.9 


2.04 



NCFC [16] 20 290.8(2) 131.5(3) 159.3(3) 2.032(2) 



TABLE VI: Results for calculations of 20 neutron drops in 
harmonic potentials with the same conventions as Table IVl 



Approx. 


ho, 


-Etot 


-Bint 


f/ext 


f rms 


HF 


5 


230.4 


120.9 


109.6 


4.26 


HF/NV 


5 


232.8 


120.2 


112.5 


4.32 


HF/PSA 


5 


231.3 


120.7 


110.6 


4.28 


BHF/PSA 


5 


221.3 


112.4 


108.9 


4.25 


fit/PSA 


5 


223.0 


117.5 


105.5 


4.18 


HF 


10 


455.4 


224.0 


231.5 


3.10 


HF/NV 


10 


466.5 


224.2 


242.2 


3.17 


HF/PSA 


10 


459.3 


224.1 


235.2 


3.12 


BHF/PSA 10 


445.0 


215.7 


229.3 


3.08 


fit/PSA 


10 


441.5 


214.8 


226.7 


3.07 


NCFC [8] 


10 


< 452. 








HF 


15 


693.0 


328.1 


364.9 


2.59 


HF/NV 


15 


715.7 


332.5 


383.2 


2.66 


HF/PSA 


15 


700.9 


329.5 


371.4 


2.62 


BHF/PSA 15 


680.1 


318.2 


361.8 


2.58 


fit/PSA 


15 


675.9 


313.2 


362.7 


2.59 


NCFC [8] 


15 678(8) 322(10) 356(10) 2.56(4) 


HF 


20 


941.3 


435.2 


506.1 


2.29 


HF/NV 


20 


976.1 


446.3 


529.8 


2.34 


HF/PSA 


20 


953.7 


438.4 


515.4 


2.31 


BHF/PSA 20 


928.1 


417.7 


510.4 


2.30 


fit/PSA 


20 


924.4 


414.5 


509.9 


2.30 



NCFC [8] 20 922(6) 425(10) 497(10) 2.27(3) 
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FIG. 1: Total energy and radius for 8 and 20 neutrons in a har- 
monic potential with oscillator parameters 10 MeV, 15 MeV 
and 20 MeV. Calculations using the NCFC are compared to 
full HF, the DME/PSA approximation to HF with exact 
Hartree, and results incorporating the resummed ladders for 
neutron matter using a density-dependent adjustment (BHF) 
and using fit coefficients from Eq. (|15[1 , 



to HF shows that the DME-HF/NV and DME-HF/PSA 
calculations have consistently higher energies and radii. 
These trends and the comparison to the NCFC results 
are evident in Figs. [T] and where energies and radii 
from the various DME prescriptions are compared to a 
full Hartrce-Fock calculation and to NCFC calculations 
for 8 and 20 neutrons. The energies have been scaled by 
the Thomas-Fermi energy trend N^Ml [3 to remove 
the dominant dependence on N and Ml. The error bars 
from the NCFC are from the extrapolations; the exact 
results for the Minnesota Hamiltonian are expected to 
lie within these error ranges. 

The DME/PSA BHF results improve on the HF and 
DME/PSA HF results systematically for both the total 
and internal energies, and for the radii. They are within 
the errors of the NCFC in many cases. The DME/PSA 
fit result with only additional volume terms fit to neu- 
tron matter, is consistently worse than the HF results, 
particularly for the internal energy (not shown). This 
failure is consistent with Ref. [l3[, where the analogous 
prescription was found to introduce unacceptably large 
over-binding in nuclei. Here, the total energies and rms 
radii are systematically too small. But once a surface 
term is added and fit to the NCFC results for the total 
energies, excellent systematics are found for these ener- 
gies and the predicted radii. The predicted internal en- 
ergies are also improved although there are much larger 
differences from the central NCFC values than for the 
total energies. 
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FIG. 2: Internal energy and radius for 8 and 20 neutrons 
in a harmonic potential with oscillator parameters 10 MeV, 
15 MeV and 20 MeV as in Fig. Q] 
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FIG. 3: Density for 8 neutrons in a harmonic potential 
with oscillator parameter 10 MeV. The non-interacting den- 
sity (HO) is compared to calculations using the NCFC 
(shaded region), the DME/PSA approximation to HF with 
exact Hartree, and results incorporating the resummed lad- 
ders for neutron matter via a density-dependent adjustment 
(DME/PSA BHF) and using fit coefficients. 



The reproduction of coordinate-space densities and the 
corresponding momentum-space form factor are shown 
for some representative cases in Figs. [3j [4j and These 
examples illustrate the very large range in interior density 
probed by this set of external potentials. The improve- 
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r [fm] 

FIG. 4: Density for 20 neutrons in a harmonic potential with 
oscillator parameter 15MeV, see Fig. [3] caption. 




q [MeV] 

FIG. 5: Form factor for 8 neutrons in a harmonic potential 
with oscillator parameter 10 MeV, see Fig. [3] caption. 

merits noted for DMA/PSA BHF energies and radii are 
also seen for the densities. The DME/PSA fit results 
might be judged the best, although the error bands from 
the NCFC calculation are too large to allow a definitive 
conclusion. The form factor is given by: 

If o , . sin qr 

F(q) = -J dr4nr 2 p(r)—^ . (16) 

Not surprisingly, the characteristic features of the first 



minimum and height of the second maximum are well 
reproduced within the NCFC error band. 



IV. SUMMARY 

In this paper, we perform test calculations to help de- 
velop density functionals for nuclei using microscopic in- 
put. In particular, we use the density matrix expansion 
(DME) as a bridge from many-body perturbation theory 
to EDF's that can be used in solvers and with optimiza- 
tion techniques developed for phenomcnological Skyrmc 
functionals. There are many implementation questions 
and options for the DME, some of which are addressed in 
the present work. Ultimately we will use high-precision 
two- and three-body nuclear interactions, such as from 
chiral effective field theory, evolved to softer forms using 
renormalization group methods. This softening makes 
them suitable for a MBPT treatment, unlike conventional 
interactions. As an interim step, we have used the Min- 
nesota potential as a (moderately) soft, semi- realistic in- 
teraction for our tests. 

The test environment is interacting neutrons in a har- 
monic trap. By varying the oscillator frequency of the 
trap, different inhomogencous density profiles are probed. 
According to DFT, an EDF for self-bound nuclei should 
be the same for the neutron systems, with the simple ad- 
dition of the external potential to the Kohn-Sham poten- 
tial and U cx t to the energy functional. Thus we can make 
controlled explorations of the energies for different den- 
sity distributions. The predictions are validated against 
ab initio calculations using NCFC methods, made pos- 
sible for 8 and 20 neutrons by computational and algo- 
rithmic advances enabled by the UNEDF project. A key 
feature is that the same Hamiltonian is used for the ab 
initio solution and the functional. This allows us to more 
reliably isolate different sources of approximation errors. 

Comparisons were first made between the DME at the 
HF level and full HF calculations. The best results for the 
improved form of the DME, which uses phase-space av- 
eraging (PSA), were consistently superior to the original 
Negele-Vautherin formulation. This agrees with the find- 
ings in Ref. [27j , where it was found in non-self-consistent 
calculations of nuclei that the PSA-DME was the most 
accurate formulation when the expansion is truncated 
at second order in gradients. (Note: the improvement 
will be more significant for more realistic interactions in- 
cluding spin-orbit and tensor forces.) Several options for 
treating the Hartrcc (direct) contribution were consid- 
ered, with a clear preference for an exact treatment in 
the PSA version of DME. This is consistent with past 
experience applying the DME. 

There are still systematic discrepancies between the 
best DME Hartree-Fock calculations and the full HF 
results. This is in contrast to results from recent 
exact-exchange orbital-based calculations, which find ex- 
tremely close agreement with HF for the same neutron 
systems and interactions [45[. Thus the source of the 
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discrepancies seems to be the DME approximation and 
not simply the use of a local Kohn-Sham potential. This 
implies that DME-based functionals will always need to 
be supplemented with some correction mechanism for in- 
herent DME approximations as well as for errors from 
many-body approximations. 

Hartree-Fock is a good quantitative starting point for 
these neutron drop systems with the Minnesota poten- 
tial, but the full NCFC solutions show clear differences 
in the patterns of energies and radii. These differences 
are used as a testbed for two ways to incorporate correla- 
tions beyond HF in the DME functional. In both cases, 
Brueckner-Hartree-Fock (BHF) calculations of neutron 
matter with the Minnesota potential were used as the 
exact reference. (Instead, one could use ab initio calcu- 
lations if available, or else calibrate to phcnomenological 
values.) 

In the first approach, the functional was modified by 
density-dependent terms, while in the second approach, 
the interaction was supplemented with Skyrme-like con- 
tact terms, whose coefficients were fit to the neutron 
matter calculation. The first method (labeled BHF in 
the figures), showed systematic improvement from DME 
Hartree-Fock toward the full NCFC results, both for the 
energies and the radius. The second method with vol- 
ume terms only was found to be unacceptable, but the 
addition of a fit surface term improved the systematics 
dramatically, to roughly the error limits of the current 
NCFC calculations for the total energies and the pre- 
dicted radii. These results validate the strategics planned 
for more realistic forces. 

There are many ways forward from the present cal- 
culations. The neutron drop system will continue to 



be a valuable tool for diagnostic testing, which will in- 
clude using non-harmonic (e.g., Woods-Saxon) and de- 
formed traps. On-going development of the density ma- 
trix expansion includes the extension to pairing, which 
will be tested using open-shell neutron numbers, and ex- 
tensions to incorporate higher-order MBPT and more 
complete low-momentum potentials. These extensions 
will be tested both in the neutron drop systems and for 
ordinary self-bound nuclei. Work in these directions is in 
progress. 
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